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ABSTRACT 

A unique formulation of describing fluid motion is presented. The method, referred 
to as “extended Lagrangian method,” is interesting from both theoretical and numerical 
points of view. The formulation offers accuracy in numerical solution by avoiding numerical 
diffusion resulting from mixing of fluxes in the Eulerian description. Meanwhile, it also 
avoids the inaccuracy incurred due to geometry and variable interpolations used by the 
previous Lagrangian methods [1,2]. Unlike the Lagrangian method proposed in [3,4] which 
is valid only for supersonic flows, the present method is general and capable of treating 
subsonic flows as well as supersonic flows. The method proposed in this paper is robust and 
stable. It automatically adapts to flow features without resorting to clustering, thereby 
maintaining rather uniform grid spacing throughout and large time step. Moreover, the 
method is shown to resolve multi-dimensional discontinuities with a high level of accuracy, 
similar to that found in one-dimensional problems. 

1. INTRODUCTION 

It is well known that fluid motion can be specified by either the Eulerian or Lagrangian 
description. Most CFD developments over the last three decades have been based on the 
Eulerian description and considerable progress has been made. In particular, the upwind 
methods, inspired and guided by the work of Gudonov[5], have met with a great deal of 
success in solving fluid flows, especially where discontinuities exist. However, this shock 
capturing property has proven to be accurate only when the discontinuity is aligned with 
one of the grid lines since most upwind methods are strictly formulated in a one-dimensional 
framework and only formally extended to multi-dimensions. Consequently, the attractive 
property of crisp resolution of these discontinuities is lost. Even though research on genuine 
multi-dimensional approaches has recently been undertaken by several leading researchers, 
they are nevertheless still based on the Eulerian description. 

Recently, Loh and Hui[3] have convincingly demonstrated that a Lagrangian formu- 
lation can capture a contact discontinuity crisply, which is difficult to achieve by Eulerian 
formulation without resorting to some special treatment such as sub-cell resolution. Fur- 
ther developments have been carried out by Loh and Liou to solve real gas problems[6] 
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and three-dimensional supersonic problems[4]. The 3D extension is not so trivial and in 
fact it involves somewhat tricky definition of cell movement and approximate solution of 
the multi- dim ensional Riemann problem. Several interesting 3D problems that have not 
been attempted previously were calculated and again as in the 2D case, high accuracy 
was achieved in resolving very complex shock-shock interactions. This method employs 
the point of view that it strictly follows the fluid particles released at some initial time 
line. The streamlines become a “time-like” coordinate and are used also for identifying 
particles. Therefore the method is naturally suitable for supersonic steady flow. No grid 
generation is needed a priori since the grid is a part of the solution, viz new grid lines are 
obtained as the solution marches in the “time-like” direction. Unfortunately, restriction to 
supersonic flows only limits the use of the method. To include the subsonic regime requires 
substantial conceptual changes. 

Numerically solving subsonic flows using the strict Lagrangian concept becomes an 
excessive obstacle. Numerous researchers at Los Alamos, making substantial contribu- 
tions to this subject, have proposed the so-called Arbitrary Lagrangian-Eulerian(ALE) 
method, for example [1,2] and others. The method consists of three phases of numerical 
procedures, using time-splitting approximations, and can become rather complicated and 
tedious. Tracking the Lagrangian cells requires interpolations of data and coordinates. 
Additional constraints must be imposed to prevent grid lines from crossing each other 
during the rezoning. Loss of accuracy is then inflicted through the continuous geometry 
and flow interpolations in response to the fluid particle motion. Indeed, such a Lagrangian 
procedure finds some analogy in the shock fitting procedure, although the former is more 
complicated. 

Physically fluid particles seem to adjust to motion and to surrounding (geomet- 
ric or physical) constraints quickly and graciously, in particular sensing the upstream- 
propagating influences. Thus, a key to the success of a numerical Lagrangian procedure 
lies in how to properly and instantaneously feed these upstream-propagating waves to the 
particles, while tracking them. It is indeed a very challenging research topic that motivates 
us to begin this exploratory investigation. This paper presents the salient features of the 
method, referred to as “extended Lagrangian method”. For flows at all speed regimes 
including purely subsonic and mixed flows, we demonstrate the advantages of the method 
over the Eulerian description, with focuses on important features commonly seen in com- 
pressible flows, such as shocks, expansion waves, slip surfaces, and interactions among 
them. We list the distinct features resulting from the extended Lagrangian method. 

1. The solution adapts to the flow variations (smooth or sudden), notably shocks and 
contacts, and as such it can be regarded as an automatic procedure for solution adap- 
tation. 

2. Unlike the conventional adaptation techniques, there is no need for clustering grid 
points near the discontinuities. Very uniform grid can be maintained and in fact can 
also achieve orthogonality easily by construction. Thus discretization accuracy does 
not deteriorate. Since streamlines do not cross, grid singularity or negative volume 
certainly will not occur. 
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3. As will be seen later, the shock capturing quality in 2D is comparable to that found 
in the ID problem. This suggests that the present approach can be viewed as an 
alternative to the current genuine multi-dimensional approach. 

4. The contact discontinuity can be resolved crisply, since it is a streamline and as such 
no numerical diffusion is introduced due to fluid crossing the cell face. 

The rest of the paper is organized as follows. In Section 2 we compare differences of 
Eulerian and Lagrangian descriptions of fluid motion, with an emphasis on the numerical 
aspects. Some current Lagrangian approaches are commented on in Section 3. Section 
4 outlines the key elements in formulating the present extended Lagrangian method for 
solving subsonic flows by retaining the advantageous features of the Lagrangian approach. 
A detailed formulation is then given in Section 5, including the discretization method 
and boundary conditions. Section 6 describes the grid motion of the present “extended 
Lagrangian” method. Finally, in Section 7 we demonstrate the advantages of the proposed 
method by displaying solutions of flows at all speed regimes, containing various interesting 
features. 


2. EULERIAN VS LAGRANGIAN DESCRIPTION 

By definition, the Eulerian description oberves at fixed locations the flow properties as 
fluid particles pass by. This has a close relation to the conventional computation approach 
in that each fixed grid point can be thought of as an observing station — corresponding to 
probes in measurement. With this approach, the meshes are generated mainly based on 
the geometry constraint, with little regard given to the motion of fluid. Naturally, the grid 
lines will seldom coincide with fluid path lines. Even when grid lines are clustered near 
high-gradient regions using conventional adaptation techniques, they are not aligned with 
the particle path. A good example is the shock wave along which grid lines are densely 
distributed, and with which streamlines make a nonzero angle, since the fluid will always 
pass through, not along, the shock wave. The angle is usually oblique in multidimensional 
flows. Consequently, the Eulerian approach has several severe numerical effects on the 
solution accuracy: 

1. Fluid particles are free to cross the grid line, thereby bringing (convecting) with them 
numerical mixing and diffusion across the cell interface. 

2. This numerical diffusion is only associated with the error resulting from approximating 
the convective terms. 

3. A contact /slip or shear layer is smeared ever increasingly with time and distance, 
unless some detection and special treatment techniques are employed. See [7] for 
example. 

In spite of numerical diffusion resulting from approximation of convective terms, the 
Eulerian description does offer convenience and simplicity both conceptually and geomet- 
rically. The grid can be constructed regularly, simply based on geometry constraints and 
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independent of flow features. To enhance accuracy, grid adaptation is often applied to 
regions of high gradients. However, the concept of adapting to high gradients inevitably 
results in a skewed and distorted grid. This feature will become more troublesome as two 
or more high- gradient areas intersect. 

It is quite safe to say that during the last three decades CFD algorithm researchers 
have primarily concentrated on developing better (more accurate/robust /efficient) ways 
to deal with the convective terms, which exist only in the Eulerian formulation. Conse- 
quently, much success has been achieved, perhaps to the point of near perfection and little 
return could be gained. Unfortunately, inaccuracy due to numerical diffusion(mixing) in 
forming the interface numerical fluxes still exists and becomes more exaggerated in multi- 
dimensional problems. On the other hand, since the convective terms do not explicitly 
appear in the Lagrangian formulation, the numerical mixing automatically disappears in 
the flux evaluation, rendering the Lagrangian approach more attractive. However, some 
other technical barriers have surfaced and discouraged researchers from further pursuing 
development of methods for this approach. In what follows, we shall detail the concept 
of the Lagrangian approach from the viewpoint of numerical solution. The differences 
between the two descriptions will then become obvious. 

The Lagrangian description, by definition, states the motion and properties of given 
fluid particles as they travel to different locations. In particular, since the particle path 
in steady flow coincides with the streamline, no fluid particles will cross the streamline. 
In other words, while staying in contact, neighboring streams will not mix via convection, 
except in the molecular level where the physical molecular diffusion takes place. Therefore, 
the following numerical consequences can be realized: 

1. No numerical diffusion is introduced across the cell interface since the computational 
cells follow the streamline. 

2. Fluid particles change motion (direction/speed) only as warranted, e.g., as shock or 
exp ansi on waves are encountered. In other words, streamlines will bend, converge, or 
diverge only as situations demand. 

3. This description gives a realistic depiction of flow behavior; cells of same j— index 
form a streamline that is identifiable with flow visualization. 

The notion of using a Lagrangian approach to describe flow is not new. In fact, 
the very essence of following fixed, particles also presents mathematical complexity to the 
approach, thereby limiting its scope of success. With the help of the new Lagrangian 
formulations, numerical solution can be as easily obtained as for the Eulerian approach, 
only with the additional distinct advantages as stated above. In the following, we shall first 
review some current numerical procedures based on the Lagrangian approach, commenting 
about their strengths and weaknesses. Then we will focus on the applicability to the more 
difficult problem, namely the subsonic regime. 
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3. REVIEW OF CURRENT L AGRAN GIAN APPROACHES 


The Arbitrary Lagrangian-Eulerian Technique (ALE) perhaps is the most well-known 
Lagrangian formulation in use at present time. The technique, initially conceived and 
developed at Los Alamos during the 70s and further implemented in the production codes 
such as CAVEAT and KIVA[8,9], etc., consists of three phases of numerical procedures 
using the time splitting concept. For a complete description of the procedure, the reader is 
referred to Refs. [1,2, 8, 9]. Continuous rezoning is carried out in order to follow the particle 
motion. As a result, spurious error produced by this procedure can lead to grid irregularity 
and tangling. Thus loss of accuracy, manifested as numerical diffusion, is inflicted through 
geometry and flow variable interpolations. 

Recently a new Lagrangian Formulation was proposed by Hui and Van Roessel [10] 
The invisicd conservation laws are transformed by using stream functions and Lagrangian 
time as independent variables. The stream functions serve to identify particles, while La- 
grangian time represents time-like coordinate. Under this formulation, geometry conserva- 
tion is enforced and each cell is literally a fluid particle. Since there is no need for remap- 
ping, the associated loss of accuracy seen in the ALE method does not appear, allowing 
extremely sharp resolution of contact discontinuities. Successful demonstrations have been 
made by Loh and Hui [3] in 2D and Loh and Liou [4] in 3D problems. Multi-dimensional 
discontinuities are resolved with the same level of accuracy as their one-dimensional coun- 
terparts, indicating that the Lagrangian formulation inherently includes multi-dimensional 
flow characteristics. However, a severe limitation restricts the validity of the formulation 
[3,4] to only supersonic flows because the formulation is based on the use of the time-like 
coordinate. Thus, extention to subsonic flows based on the same framework appears im- 
possible. In what follows we will first give the basic ideas for extension in the next section 
and then describe detailed steps in Section 5. 


4. EXTENSION TO SUBSONIC FLOWS 


A key element in the subsonic flow is the existence of the upstream-propagating wave. 
Thus, the existence of a body located downstream is transmitted to the oncoming fluid 
particles via this wave so that the particles can change motion accordingly. This immedi- 
ately implies that we must abandon the time-like formulation since it is only suited for pure 
initial value problems, such as supersonic flow where no influence comes from downstream. 
Next, we must also abandon the idea of following a fixed particle, at least for the steady 
flows. Alternatively, we consider the steady streamlines as a set of lines that are occupied 
by particles released at the same location, different times and yet treated indistinguishably. 
The upstream-propagating influence is felt through the downstream particles on the same 
str eamli ne in order to satisfy the governing conservation equations and boundary condi- 
tions in question. By describing fluid motion along streamlines, we allow fluids to maintain 
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their identity without tracking each specified particle, This definition is of course broader 
than and is a sufficient condition to the Lagrangian description, which follows motion of 
fluids of spcific identity. Consequently, the present method is termed extended Lagrangian 
method. The net result is that we retain the essential beauty of the Lagrangian description 
that introduces no or minimal numerical diffusion across streamlines. 


5. APPROACH 


To facilitate the description, let us first define the notation for the relevant variables 
in the 3D Euler equations. The physical variables in a phase space of dimension 5 are 
denoted by a boldface uppercase letter or column vector whose elements are denoted by 
lowercase letters. 


U = 


f p\ 


f P \ 

pu 


pu 

pv 

, u c = 

pv 

pw 


pw 

\pe t J 


\ph t / 


( 1 ) 


where e% = e + 0.5(u^ + + w^') and ht — e* + p/ p. The geometrical vectors in physical 

(Cartesian) space of dimension 3 are denoted by an overhead arrow The fluid velocity 
is _ _ _ 

V' = ui + vj + wk, (2a) 


and the normal vector of the boundary surface of a control volume 

S = s x i + Syj + s z k. 

The inviscid fluxes in 3D physical space are compactly written as 


F = 


(26) 


l p \ 


( °-\ 


o\ 

pu 

v + 

Pi 

= U c y + P, where P = 

Pi 

g 9 
Q. CL 

°iy2. 
^ 

PI 

pk 

0 / 


( 3 ) 


The first term in F is the flux of U c convected by the fluid velocity V and the second 
term simply the pressure flux. 

For the following discussion, it is useful to review some basic concepts used to describe 
fluids. It is understood that the fluid has been considered to be a continuum. A convenient 
concept within continuum mechanics for describing a fluid motion is that of control volume. 
In Fig. 1, let fi*(/) be a moving volume with bounding surface d£l*(t)\ the local boundary 
velocity is V&. The volume is arbitrary and in general need not be identified with either 
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physical boundary or specific motion of the fluid in ft*. Such a volume is called control 
volume. A special type of control volume is called material volume, denoted by ft(t), 
consisting of a collection of matter of fixed identity enclosed by a material surface, denoted 
by 9ft(<), of which every point moves with the local fluid velocity V. (See also Fig. 1.) If 
the volume ft(t) is shrunk to a point, the resulting material volume is called a fluid particle. 
Consequently, the fluid properties of the fluid particle can be described mathematically in 
terms of space and time. 

Under the assumption of the continuum mechanics, let x be any continuous function, 
such as the density. The Reynolds transport theorem [11] gives the time rate of change of 
the “content” x of ft*: 


J xdv= J ^dv + J xVfdS, (4) 

o*(t) n*(t) an »(t) 


where the element surface dS of ft* is moving with the velocity Vj. Note that Vj may vary 
over the surface 9ft*. Again, a special case is when the theorem is applied to the material 
volume ft(<) with Vb = V. 

The conservation laws (neglecting viscous fluxes for simplicity, without loss of gener- 
ality in decribing the approach) can be conveniently expressed over an arbitrary control 
volume ft*(<) in integral form: 


j t J Vdv+ J [\J c (V-V h ) + P b ]*dS = 0, 

n-(o an-(t) 


with Pfc = P + 


/ 0 \ 

0 

0 

0 

\pV b J 


( 5 ) 


From the above equations, three fundamentally different approaches can result. 


( 5 . 1 ) Eulerian Description-. 

The Eulerian description assumes that the observer stays stationary with respect to 
the chosen frame of reference (e.g., inertial system). This requires: 

V b = 0, and ft* ^ ft*(t). (6) 

That is, the control volume is fixed in time. 

With the application of the Reynolds transport theorem Eq. (4), Eq. (5) is reverted 
to the familiar integral form: 

Jw dv+ J [VcV + P]*d 3 = 0 , (7) 
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where the superscript “E” is used to denote Eulerian frame of reference. In the discrete 
version of the above equation, each cell represents a control volume and is not moving in 
time, even though the flow may be unsteady. Note that the case in which an observer 
is fixed to a non-inertial frame, e.g. on rotating machines, also belongs to the Eulerian 
description. 

(5.2) Lagrangian Description: 

According to the strict definition, the Lagrangian description requires that the volume 
fl *(t) (enclosed by d£l*(t)) move with the instantaneous fluid velocity and is identified as 
the material volume fl(t) (see also Fig. 1). That is, 

V* > 0, and n*(<) = ft(<). (8) 

And the conservation laws are simplified to 

£ J U dv+ j P h *d§ = 0. (9) 

n(t) an(t) 

Clearly, the pressure remains as the only contribution to the flux on the surface. This 
makes it extremely simple to calculate the time-rate of change of f Udv, i.e., involving 
only the pressure acting on the bounding surfaces. However, the trajectory of the vertices 
is the part that often causes difficulties, resulting in large deformation or irregularity. 
(See [1,2].) Nevertheless, this is a very intriguing idea that avoids the nonlinearity in the 
equation of motion, thus reducing many difficulties associated with the convective terms 
that exist only in the Eulerian viewpoint. Such nice properties unfortunately have not 
been able to outweigh the drawbacks (see Introduction) and gain favorable reception vs 
the Eulerian approach. 

In the following, we propose a new approach that retains essential advantages of the 
above two approaches. 

(5.3) Extended Lagrangian Description : 

Close investigation of the surface integral reveals that the convective term can be 
eliminated also by requiring that: 


(V-V b )»S = 0. (10) 

That is, a portion of the surface is parallel to (V — V b ). Equation (5) then becomes, 

j t J Vdv+ J p b *ds+ J [(y-H)u c + p*W5 = o. (ii) 

3nf L (t) dn$ L (t) 
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The control volume now is denoted by superscript “EL” to indicate the present description. 
The surface d£l EL is comprised of two types, dQ, EL = dCl EL UdO, EL , where dQ, EL coincides 
with the instantaneous particle paths and dCl EL represents the inflow and outflow faces. 
Clearly, the present “extented Lagrangian method” combines the quality of the above two 
approaches: the second surface integral includes both convective and pressure terms, as in 
the case of Eulerian approach; the first surface integral, on the other hand, merely has the 
effect of pressure, as in the case of Lagrangian approach. That is, see Fig. 2, 

On dU EL +— both convected and pressure fluxes 

On dCl EL <— only pressure flux 

For steady flow, Vb = 0, there is no need for literally following particles because 
no variations of motion with time appear among the particles on the same streamline. 
Therefore that whether we strictly follow particles of same identity or not is irrelevant in 
the formulation. Indeed, following the streamlines, rather than particles, is the essence 
of the present approach and this rescues us from facing the difficulty of other Lagrangian 
approaches. Substituting Vb = 0, in Eq. (11), we get 


j t J U dv+ j P»dS+ J [U c y + P] • dS = 0, m9, EL . 

QEL QftEL d^EL 


(12a) 




together with the constraint, 

V’ • S = 0, on dCl EL . (12&) 

For steady flow (fixed volume), the semi-discrete form, with the time-dependent term 
retained for iteration purpose, can be cast as: 


- f 

dt J 


Vdv + 53 ft # 4 + 53 ft •4 = 0. 

iednf L iea n* £ 


(13) 


Ex ami nation of the above equations reveals some interesting insight. The inbalance of 
pressure in two neighboring cells with common interface boundary dO, EL causes the change 
of flow direction(i.e., dQ, EL ) of the fluids under consideration sis well sis change of their 
volumes. In other words, the deformation and dilatation of the fluid can be described. 
Indeed, the Lagrangian grid includes multi- dimensional information and suggests how the 
fluid volume distorts in the flow. This point of view malces the description of fluid motion 
intuitively simple and clear. Moreover, it results in a major benefit in the numericsil 
solution because this formulation avoids any arbitrary(numerical) mixing of fluids which in 
turn introduces numericsil diffusion in the solution, notably across the contact discontinuity 
or shear layer. This diffusion error is common in the Eulerian formulation in which a 
nonstationary contact discontinuity is smeared without bound as time/space is marched. 
Furthermore, the advantage of the present approach is also clearly shown in its capability 
for crisply resolving multi-dimensional shocks. 
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To complete the numerical solution procedure in the finite voulume, we employ an 
accurate and efficient new upwind scheme described in full detail in [12,13]. In what 
follows, we shall see that this new splitting has a very interesting bearing with the present 
extended Lagrangian method. 


(5.4) The Upwind Method[l2,13] 

To illustrate the concept, it is sufficient to consider only the one- dimensional system. 
As a first step, by recognizing convection and pressure as two physically distinct (but 
coupled) processes, we split the fiux in the form of Eq. (3). In other words, these two 
terms deserve separate treatments. Mathematically, we propose to separately deal with 
the genuinely nonlinear ((it — a, u + a) pair) and linearly degenerate (u) fields. 


F = u | pu ] -f P = F c + P, F c = uU c . (14) 

\phtj 

The overhead arrow has been dropped for we are concerned only with one-dimensional 
flow. Both Mach n umb er and velocity splittings can be used to represent the convective 
quantity u in F c . In most cases, there is virtually no difference between calculated results 
of the two splittings. As found in [13], the velocity splitting is more robust in calculating 
unsteady shock tube problems by allowing a larger time step at start. Now, the numerical 
convective flux at the interface (denoted by subscript |) straddling the left (I) and right(i?) 
states, is effectively written as: 

F cx /2 = U 1 / 2 U cL/H, (1®) 

where Ui/ 2 is the interface convective velocity. Let Ui/ 2 be written as: 

U1/2 ~ Ur- (I®) 

Several formulas are appropriate to define u ± , e.g., 


v ± = J ( u ± M)/ 2 > 

\ ±(tt ± a) 2 /4a, 


if |u| > a, 
otherwise, 


(17) 


where a is the speed of sound. The convectible variable vector U c is then upwinded solely 
based on the sign of u x / 2 , viz, 


TT - / ( U c)^ if U 'l* * °’ 

cL/H ^ (u e )*, otherwise, 


(18) 


We turn now to the pressure term by writing: 

P1/2 — Pl +Pr- 


(19) 
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Similarly, a whole host of choices are possible for the pressure splitting. A differentiable 
pair of the ‘+’ and components have been found to be effective. 

„± _ / Pi 1 ± sgn(u))/2, if |u| > a, /gO') 

** \p(M ± 1) 2 (2 ^ M)/4, otherwise. ' 

This completes the definition of the numerical flux F. Putting (15) and (19) together, 
we recast the interface flux in the following form: 

F 1/2 = Ul/2 ~(U C £ ^ 2 l Ul / 2 l ^l/2^c + P 1/2' (21) 

where Ai/ 2 {*} = {*}r — {•} Here the first term on the RHS is clearly not a simple 
average of the ‘L’ and l R > fluxes, but rather a weighted average via the convective velocity. 
The dissipation term has merely a scalar coefficient |iti./ 2 1 and requires only 0(n) oper- 
ations for an n— dimensional vector F. Furthermore, since there is no differentiation (or 
jacobian matrix) involved in evaluating F^, the present method is easily extended to a 
general equation of state and non-equilibrium flows and the cost is only linearly increased 
with the additional conservation equations considered. Unlike the Roe or Osher schemes, 
the extension does not yield additional ambiguity such as the definition of averaged or 
intermediate states. Also, numerical tests strongly suggest an entropy-satisfying property 
by the present method. 

To achieve higher-order spatial accuracy, a MUSCL-type procedure is followed to 
upwind-extrapolate variables(primitve variables in the calculations presented in this pa- 
per), with TVD limiters incorporated [14]. Then, a two-stage Runge-Kutta procedure is 
used to integrate the semi-discrete sytem Eq. (13), subject to the kinematic condition Eq. 
(12b). 

The subsonic inflow conditions are imposed by specifying total enthalpy, total pressure, 
and flow angle, while the outflow conditions are obtained with specified static pressure and 
extrapolated total enthalpy, total pressure, and flow angle. The usual tangency procedure 
is used at the cell boundary that coincides with a physical wall — no ghost cells are used. 
The wall pressure is gotten using linear extrapolation from interior data, so is the total 
enthalpy. 


6. MOTION OF LAGRANGIAN GRIDS 


An important integral part of the present method is the grid motion that follows 
the constraint Eq. (12b). Two basic settings can be chosen for defining the motion of 
computational cells, namely the motion of cell centers or cell vertices. With the former 
approach the cell vertices will be defined by the position of neighboring centers, vice versa 
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for the latter. Since the constraint, Eq. (12b), is imposed on the cell boundary, it is 
consistent to determine the vertex motion instead. This is easily done with the velocity 
field known from the solution. The constraint Eq. (12b) is equivalent to the kinematic 
condition on a streamline: 

dx = dy = d 1 (22) 

u V w 

As flow variables are defined at cell centers, the velocity components at cell boundary must 
be defined by some interpolation procedure from surrounding cells. In the present report, 
we outline the general notion of grid movement and give a specific strategy showing how the 
grid is moved to meet the constraint for the test cases included in the paper. Let us consider 
the two-dimensional cell shown in Fig. 3. Three dimensional cells can be treated 

similarly. Assuming the cell boundary is described by a fine segment (rj+i j-t-i — r, -j.fi). 
Since the segment is a part of a streamline, Eq. (22) gives 

y«+i,i+i = Vi,j + 1 + m(r l+u+1 - *ij+i) (23) 


where 


m = 


v 

u 


Here we list some possibilities for evaluating (u, v) : 


(a) Mid-point average 


V = |[2$j + 2Kj+i + Vi+ij + K+ij+i + Vi-ij + K_ij+i], 


(24a) 


(b) Upstream average, assuming Ufj > 0, Vi, j, 

V = -[Kj + Kj+i + K-ij + K-_ij.fi]. (246) 

There are two unknowns, (x,y),+i j+i, in Eq. (23). Another condition is needed to 
complete the system. In this report, we prescribe the value of x-coordinate for each ith grid 
line, i.e., x =constant lines. This condition provides simplicity, but also yields accuracy as 
will be shown later. Furthermore, specification of x-coordinate allows one to put fine grids 
to resolve geometry details, e.g., near high curvature region. 

When the constraint, Eq. (12b), is satisfied for all jth grid line, the conservation laws 
are basically solved in a one-dimensional stream tube, because there is no flow across the 
j— grid lines. As a result, high accuracy is expected with this virtually one-dimensional 
problem. This is the reason that the present method gives sharp representation for oblique 
shocks, as accurate as their counterpart in one-dimension. The formulation itself already 
inherits multi-dimensional information via the deformation of grid lines (i.e., streamlines). 
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The cost of arriving at the above constraint is negligible even if it is done at each 
iteration, because calculation of Eqs. (23) and (24) is all it needs. The grid motion can be 
predicted in phase as the evolution of the flow variables, or else for a prescribed number of 
flow variables iterations. However, it is unnecessary to adjust the grid so frequently while 
the flow is still evolving. A more thorough investigation about the optimal iteration per 
Lagrangian grid motion and its sensitivity to flow condition is useful, but is beyond the 
scope of the present paper. 


7. TEST PROBLEMS AND DISCUSSION 

We will show examples for flows at all speed regimes, featuring solution accuracy and 
grid aspects. The plots are organized uniformly for all cases presented. We show Mach 
contours overlaid on the grid used. Fine grids, by doubling grid number in both direction, 
are used only in the Eulerian calculation for comparison purpose. The grids shown are 
only one half of the whole grid in the Lagrangian case, and one quarter in the fine-grid 
case, thus corresponding to roughly the same location in both plots. The symbols denote 
the entended Lagrangian solutions and the lines are the Eulerian solutions. The Mach 
contours are chosen for presentation so that any numerical anomaly or inaccuracy can be 
more easily depicted than the other variables such as pressure. 

The first example is the purely subsonic flow in which an = 0.4 flow enters a 
channal with 20% circular bump, see Fig. 5 for detailed geometry. The Mach contours, 
given in Fig. 4, show nearly perfect symmetry about the midchord, except in the wake 
region. The- wake region suggests an entropy production( numerical diffusion), likely due 
to numerical wall boundary condition, which to my knowledge is still a gray area in CFD. 
While this may be the situation, the Lagrangian solution still definitely results in a narrower 
wake, roughly half the width of the Eulerian result. The Eulerian solutions in fact have a 
slight asymmetry near the top wall, even though the residuals were dropped to machine 
zero. Notice that the computation domain is considered to be small, extending one chord 
length upstream and downstream from the bump, for this purely subsonic problem. It is 
worth noting that the present grid automatically evolves from initial Eulerian grid into the 
grid system seen in Fig. 4, according to Eqs. (23) and (24a). The grid spacing between 
streamlines increases near the stagnation points and converges as flow accelerates. The 
detailed distribution of variables on the top and bottom walls are plotted in Fig. 5. The 
fine-grid Eulerian solutions are included for comparison. The agreement is remarkable and 
again symmetry is quite evident. However, the fine-grid Eulerian solution over- predicts the 
stagnation pressure, whose theoretical value is 1.116, while the Lagrangian solution closely 
matches. Hereafter, for briefness and contrast to the Eulerian solution, we shall take the 
liberty of loosely using the term “Lagrangian solution” to mean the solution obtained by 
the extended Lagrangian method described in this paper — not in the strict Lagrangian 
sense. It is also noted that the fine- grid solution took considerably more iterations than 
the Lagrangian solution to converge. 
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The second example involves the popular test of transonic flow in a channel with 10% 
hump and Mqo = 0.675. Results axe given in Figs. 6 and 7. Again, we show the Mach 
contours and distributions on both walls. The agreement of the coaxse-grid Lagrangian 
solutions with the fine-grid Eulerian solutions is excellent. The shock resolution from the 
Lagrangian method is outstanding, so is the prediction of the so-called w Zierep” singularity 
at the foot of the shock on the curved surface. The Lagrangian solution again yields a grid 
that senses the global flow characteristic. Notice that no clustering of grid is necessary to 
capture the shock in the correct location. 

The next case is also a standard one, involving an Moo = 1.4 flow and 4% bump. 
This case consists of small subsonic pocket resulting from short Mach stems on both walls, 
shown in Fig. 8. Shock- shock-expansion waves interactions take place behind the trailing 
edge. The shock locations axe in excellent agreement with the fine-grid solutions. Close 
examination of the Lagrangian grid shows that the gnd lines (also streamlines) remain 
straight until the shock is encountered, as it should. The Eulerian grid however already 
began bending at the leading edge for all j-th lines, simply because of geometry constraint. 
The present method yields grids that sire conforming with the flow features by bending, 
expanding / contract ing . The net result is that excellent shock resolution is obtained, even 
though the shock is oblique to the grid line. Figure 9 displays a one-cell capturing of the 
oblique shock. This is not entirely surprising since the present formulation has already 
taken account of the multi-dimensional nature of the flow via the streamline deformation 
caused by the fluxes (pressure forces) of surrounding fluids. 

The fourth test is an Moo — 1-8 flow over a 15° ramp. This case consists of a Mach 
stem about 10% of channel height, a slip line emanating from the triple point, and reflected 
shocks. In Fig. 10, the mach contours depict an overall picture of the flow, demonstrating 
a sharp resolution of the ramp shock, Mach stem, and the subsequent shocks. The slip line, 
whose strength is being weakened by the expansion wave generated at the ramp shoulder 
and transmitted through the first reflected shock, is resolved to the same level of accuracy 
as given by the fine-grid solution, i.e., with the same level of spatial spreading. Figure 11 
vividly displays how the grid lines bend as the shock is encountered and change direction 
according to the flow. The grid itself already suggests the flow structures, train of shock 
reflections, expansions, as well as the Mach stem across which there is no change of flow 
angle. It is worth noting the clear slipline emanating from the triple point. In contrast to 
the shock-aligned grid, the present grid is aligned with the streamlines, which will never 
cross each other, but the shocks can. Thus the present method is indifferent to whether the 
high-gradient regions intersect. The profiles (Fig. 12) on the walls show good agreement 
of both solutions. Excellent shock resolution capability is observed on both walls even the 
second reflected shocks remain well resolved. 

CONCLUDING REMARKS 

We have presented a unique formulation for dealing with subsonic as well as supersonic 
flows. The method, referred to as extended Lagrangian method, combines the accuracy 
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belonging to a Lagrangian description and the robustness and simplicity of an Eulerian 
description. Through systematic comparison with fine-grid solutions and a theoretical 
check for flows at various regimes, we have demonstrated its capability for crisply capturing 
high-gradient regions that are not aligned with the grid. In contrast to adaptive approaches 
reported to date, the present approach already automatically inherits the ability to adapt 
to flow features, but based on an entirely different adaptive philosophy in which there is no 
need for clustering grid lines. Since the grid lines are not aligned with high-gradient areas, 
they are maintained regular and uniform. Moreover, one set of grid lines that coincides 
with stre amli nes depicts vividly a form of “numerical flow visualization”. Without resorting 
to arbitrary detecting criteria, the present approach not only predicts well the nonlinear 
waves, such as shock and rarefaction waves, but also is especially amenable to treating 
a lin early degenerate field, such as a contact discontinuity. Furthermore, since the grid 
spacing is maintained relatively unform, a large time step is permitted throughout the 
calculation, thus increasing efficiency. Also the common adaptive strategy will have great 
difficulty in the case of intersecting shocks, because the grid lines will cross each other, if 
not checked. We also suggest that the present extended Lagrangian method is a viable 
alternative approach to the current multi-dimensional scheme and grid-enriching adaptive 
procedure for complex flows having high-gradient regions. In fact it is an elegant and 
effortless approach to deal with multi-dimensional flows. Further development and 3D 
applications are currently underway and will be reported in the future. 


REFERENCES 

1. Hirt, C. W., Amsden, A. A., and Cook, J. L., “An Arbitrary Lagrangian-Eulerian 
Computing Method for all Flow Speeds,” J. Comput. Phys. 14, 227 (1974). 

2. Pracht, W., “Calculating Three-Dimensional Fluid Flows at All Speeds with an Eulerian-| 
Lagrangian Computing Mesh,” J. Comput Phys. 17, 132 (1975). 

3. Loh, C. Y. and Hui, W. H., “A New Lagrangian Method for Steady Supersonic Flow 
Computation I. Gudonov Scheme,” J. Comput Phys. 89, 207 (1990). 

4. Loh, C. Y. and Liou, M.-S., “A New Lagrangian Method for Three-dimensional Steady 
Supersonic Flow,” submitted for publication. 

5. Gudonov, S. K., “A Difference Method for the Numerical Calculation of Discontinuous 
Solutions of Hydrodynamic Equations,” Mat Sb. 47, 271 (1959). 

6. Loh, C. Y. and Liou, M.-S., “Lagrangian Solution of Supersonic Real Gas Flow,” J. 
Comput. Phys. 104, 150 (1993). 

7. Harten A., “ENO Schemes with Subcell Resolution,” J. Comput Phys. 83, 148 
(1989). 

8. Addessio, F. L. et al, “CAVEATtA Computer Code for Fluid Dynamics Problems with 
Large Distortion and Interned Slip,” Los Alamos National Laboratory, LA- 10613-MS, 
1986. 


15 



9. Amsden, A. A., Rourke, P. J., and Butler, T. D., “KIVA-II:A Computer Program 
for Ch emi cally Reactive Flows with Sprays,” Los Alamos National Laboratory, LA- 
1 1560-MS, 1989. 

10. Hui, W. H. and Van Roessel, H., “Unsteady Three-Dimensional Flow Theory via Ma- 
terial Functions,” in AGARD Symposium on Unsteady Aerodynamics-Fundamentals 
and Application to Aircraft Dynamics, Gotting, West Germany, 1985, SI, CP-386. 

11. Thompson, P. A., Compressible-Fluid Dynamics, p.15, McGraw-Hill, 1972. 

12. Liou, M.-S. and Steffen, C. J., “A New Flux Splitting Scheme,” NASA TM 104404, 
1991; also to appear in J. Comput. Pkys. 

13. Liou, M.-S., “On a New Class of Flux Splittings,” Lecture Notes in Physics 414, 
pp.115-119, Springer- Verlag, 1993. 

14. Liou, M.-S., and Hsu, A. T. “A Time Accurate Finite Volume High Resolution Scheme 
for Three Dimensional Navier-Stokes Equations,” AIAA paper 89-1994-CP, 1989. 


16 



dti(t + At) 



dtt(t) 


Control volume, denoted by may be Material volume, V b = V 
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Fig. 1. D efini tion of a control volume and material volume. 
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4. Maxi contours for — 0.4, 20% bump flow; (a) Lagrangian 
solution (134 x 70 grid), (b) Eulerian solution (266 x 138 grid). 
(Contour levels: M m ,- n = 0.0, M max = 0.8, AM = 0.02.) 


19 


p/p 



20% BUMP ON A CHANNEL 
= 0.4 

0(h z ), 134x70 & 266x1 38( ) 

a Top wall, ° Bottom wall 




Fi* 5 Distributions of variables on top and bottom walls for - 0.4, 20% bump 
' flow; lines are Eulerian solution (266 x 138 grid) and symbols are Lagrangian 
solution (134 x 70 grid) with A and O denoting top and bottom walls respec- 

tively. 20 



Fig. 6. Mach contours for M ^ = 0.675, 10% bump flow; (a) Lagrangian 
solution (134 x 70 grid), (b) Eulerian solution (266 x 138 grid). 
(Contour levels: M min = 0.0, M max = 1.6, AM = 0.04.) 
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Fig. 7. Distributions of variables on top and bottom walls for = 0.675, 10% bump 
flow; lines are Eulerian solution (266 x 138 grid) and symbols axe Lagrangian 
solution (134 x 70 grid) with A and O denoting top and bottom walls respec- 
tively. 
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Fig. 8. Match contours for = 1.4, 4% bump flow; (a) Lagrangian 
solution (134 x 70 grid), (b) Eulerian solution (266 x 138 grid). 
(Contour levels: M min = 0.4, M max = 2.2, AM = 0.045.) 
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Fig. 9. Distributions of variables on top and bottom walls for M 0 0 = 1.4, 4% bump 
flow; lines axe Eulerian solution (266 x 138 grid) and symbols are Lagrangian 
solution (134 X 70 grid) with A and O denoting top and bottom walls respec- 
tively. 24 




Fig. 10. Mach contours for = 1.8, 15° ramp flow; (a) Lagrangian 
solution (100 x 41 grid), (b) Eulerian solution (198 x 80 grid). 
(Contour levels: M min = 0.4, M max — 2.2, AM — 0.045.) 
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Fig. 11. Lagrangian grid, showing streamlines and the trains of shock and 
expansion waves. 
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Fig. 12. Distributions of variables on top and bottom walls for — 1.8, 15° ramp 
flow; lines are Eulerian solution (198 x 80 grid) and symbols are Lagrangian 
solution (100 x 41 grid) with A and O denoting top and bottom walls respec- 
tively. 9.R 
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